Usenet 1993 July/InfoMagic USENET CD-ROM July 1993.ISO/sources/unix/volume25/mrandom

Text File | 1991-12-13 | 29.8 KB | 880 lines

Newsgroups: comp.sources.unix
Path: pa.dec.com!vixie
From: [e-mail] (Paul Vixie)
Subject: v25i023: A random number generator with persistent state
Message-ID: <[e-mail]>
Sender: [e-mail] (News)
Organization: DEC Palo Alto
Date: Sat, 14 Dec 91 03:34:07 GMT
Approved: [e-mail]
Lines: 868

Submitted-By: Clark Thomborson <[e-mail]>
Posting-Number: Volume 25, Issue 23
Archive-Name: mrandom

I wrote this package to overcome some troubles I had with the random()
package.  I had been saving random()'s state table to a disk file, then
restarting the random sequence in the next program run.  I discovered some
seriously non-random behavior in the numbers resulting from this practice.
Further investigation (including examination of the object code for
random() -- painful!) showed me that it is necessary to count the number of
calls to random() in order to safely restart it.  Hence the enclosed code.

[ I added all, clean, and install targets to the Makefile.  --vix ]

#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#        "End of archive 1 (of 1)."
# Contents:  MANIFEST Makefile README mrandom.3 mrandom.c mrandom.h
#   mrtest.c
# Wrapped by [e-mail] on Fri Dec 13 14:28:31 1991
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'MANIFEST' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'MANIFEST'\"
else
echo shar: Extracting \"'MANIFEST'\" \(331 characters\)
sed "s/^X//" >'MANIFEST' <<'END_OF_FILE'
X   File Name        Archive #    Description
X-----------------------------------------------------------
X MANIFEST                   1    This shipping list
X Makefile                   1
X README                     1
X mrandom.3                  1
X mrandom.c                  1
X mrandom.h                  1
X mrtest.c                   1
END_OF_FILE
if test 331 -ne `wc -c <'MANIFEST'`; then
    echo shar: \"'MANIFEST'\" unpacked with wrong size!
fi
# end of 'MANIFEST'
fi
if test -f 'Makefile' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'Makefile'\"
else
echo shar: Extracting \"'Makefile'\" \(267 characters\)
sed "s/^X//" >'Makefile' <<'END_OF_FILE'
CFLAGS= -O
X
all: mrtest
X
clean:
X    -rm *.o
X    -rm mrtest
X
install:; @-echo "just link mrandom.o into your own programs as needed."
X
mrtest: mrandom.h mrandom.o mrtest.c
X    cc $(CFLAGS) -o mrtest mrtest.c mrandom.o
X
mrandom.o: mrandom.h mrandom.c
X    cc $(CFLAGS) -c mrandom.c
END_OF_FILE
if test 267 -ne `wc -c <'Makefile'`; then
    echo shar: \"'Makefile'\" unpacked with wrong size!
fi
# end of 'Makefile'
fi
if test -f 'README' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'README'\"
else
echo shar: Extracting \"'README'\" \(1662 characters\)
sed "s/^X//" >'README' <<'END_OF_FILE'
I wrote this package to overcome some troubles I had with the random()
package.  I had been saving random()'s state table to a disk file, then
restarting the random sequence in the next program run.  I discovered some
seriously non-random behavior in the numbers resulting from this practice.
XFurther investigation (including examination of the object code for
random() -- painful!) showed me that it is necessary to count the number of
calls to random() in order to safely restart it.  Hence the enclosed code.
X
I believe random(), with my modifications, is superior in many respects
to its competition within 4.3bsd Unix, namely rand() and rand48().
Those generators are based on a multiplicative congruential scheme,
which makes them difficult to use in applications where one is generating
points uniformly distributed on the unit square.  Certainly rand() has
the problem that a high-resolution 2-d plot of the x-y points generated by
X    for (i=0; i<n; i++) {
X      x[i] = rand()/MAXLONG;
X      y[i] = rand()/MAXLONG;
X    }
will show stripes.  Using rand()%m is suicidal, since its low-order bits
are known to be cyclic.  Perhaps I'm selling rand48() short, since I haven't
examined it in detail: quite possibly its constants were chosen to minimize
the separation between the stripes in a 2-d plot.
X
To help you get started, I have also enclosed a test-driver called mrtest.
X
Please let me know if you find any defects in the sequence generated
by the default seed built into mrtest.  I picked that seed for
sentimental reasons: it was the last number generated by my old,
and now discredited, random() sequence.
X
X                        Clark Thomborson
X                        [e-mail]
END_OF_FILE
if test 1662 -ne `wc -c <'README'`; then
    echo shar: \"'README'\" unpacked with wrong size!
fi
# end of 'README'
fi
if test -f 'mrandom.3' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'mrandom.3'\"
else
echo shar: Extracting \"'mrandom.3'\" \(7553 characters\)
sed "s/^X//" >'mrandom.3' <<'END_OF_FILE'
X.TH MRANDOM 3 "11 September 1991"
X.SH NAME
mrandom, frandom, init_rng, save_rng, restart_rng, reconstruct_rng, describe_rng \- wrappers for 4.3bsd random
X.SH SYNOPSIS
X.nf
X.B long  mrandom(m)
X.B long m;
X.LP
X.B double frandom()
X.LP
X.B int init_rng(seed, filename)
X.B long seed;
X.B char *filename;
X.LP
X.B int  save_rng(filename)
X.B char *filename;
X.LP
X.B int restart_rng(filename)
X.B char *filename;
X.LP
X.B void reconstruct_rng(seed, count1, count2)
X.B long seed;
X.B long count1;
X.B long count2;
X.LP
X.B char  *describe_rng(rngid)
X.B char  rngid[RNGIDSTRLEN];
X.fi
X.IX  "mrandom function"  ""  "\fLmrandom\fP \(em generate random integer mod m"
X.IX  "frandom function"  ""  "\fLfrandom\fP \(em generate random 64-bit float"
X.IX  "init_rng function"  ""  "\fLinit_rng\fP \(em initialize random generator"
X.IX  "save_rng function"  "" "\fLsave_rng\fP \(em save random generator to file"
X.IX  "restart_rng function"  "" "\fLrestart_rng\fP \(em restart random generator from file"
X.IX  "reconstruct_rng function"  ""  "\fLreconstruct_rng\fP \(em restore random generator from (seed, count1, count2) values"
X.IX  "describe_rng function"  "" "\fLdescribe_rng\fP \(em construct short string describing generator state"
X.IX  "random number generator"  "\fLmrandom\fP"
X.IX  "random number generator"  "\fLfrandom\fP"
X.IX  "random number generator"  "\fLinit_rng\fP"
X.IX  "random number generator"  "\fLsave_rng\fP"
X.IX  "random number generator"  "\fLrestart_rng\fP"
X.IX  "random number generator"  "\fLreconstruct_rng\fP"
X.IX  "random number generator"  "\fLdescribe_rng\fP"
X.IX  "generate random numbers"  "\fLmrandom\fP"
X.IX  "generate random numbers"  "\fLfrandom\fP"
X.IX  "generate random numbers"  "\fLinit_rng\fP"
X.IX  "generate random numbers"  "\fLsave_rng\fP"
X.IX  "generate random numbers"  "\fLrestart_rng\fP"
X.IX  "generate random numbers"  "\fLreconstruct_rng\fP"
X.IX  "generate random numbers"  "\fLdescribe_rng\fP"
X.SH DESCRIPTION
X.LP
X.I mrandom(m)
generates random integers in the range 0 to m-1, using a call to the
X4.3bsd random() function.
This is a non-linear additive feedback random number generator.
The main advantage of mrandom(m) over random() is that you can save the
complete state of the random number generator to a disk file,
for later continuation of the same pseudorandom sequence.
The period of mrandom(m)'s pseudorandom output is thus guaranteed to be long.
In contrast, random()'s period can be quite short if you repeatedly
save and restore its state using setstate(), because some of random()'s
state is private.
X.LP
To use mrandom(), you must first initialize the random number generator.
The easiest way to do this is to make a call to
X.I init_rng(seed, filename).
If all goes well, this call will return the value 1,
indicating that a file has been created
describing the initial state of your generator.
Hereafter, this file will be called an RNG state file.
X.LP
X.I init_rng
will perform some error checks on the RNG state file after it is written,
returning the value 0 if a problem is detected.
X.LP
After an
X.I init_rng(seed, filename)
call, the random number generator is in the same state as it would
be if you had directly called
X.I srandom(seed)
or
X.I initstate(seed, state, 128)
on a ``fresh'' copy of the random() code.
X.LP
Warning: do not write your own calls to
X.I srandom(), initstate(), random(),
or
X.I setstate()
if you are using the mrandom package.
If you do so, the long-period advantages of mrandom will be lost,
and a later call to save_rng or restore_rng may fail.
X.LP
It is good experimental practice to preserve copies of
your RNG state files, in case you ever want to reproduce a result.
The state file is ASCII text, so it is portable across systems.
X.LP
Another good experimental practice is to ``stamp'' each printout
with a record of the initial RNG state.
You can do this by including a copy of the RNG state file, although
this is a dozen lines long and a few hundred characters.
A more convenient method is to use a call to
X.I describe_rng(rngid).
This puts a one-line short-form of the current RNG state in the
character string rngid.
X.LP
X.I describe_rng(rngid)
returns the value of rngid, for convenience.
To print this short-form RNG descriptor, for example,
you can execute printf(describe_rng(rngid)).
Note: you must declare storage for rngid in your own program with
a declaration such as
X.I char rngid[RNGIDSTRLEN].
X.I RNGIDSTRLEN
is defined with the value 80 in the current release of mrandom.
X.LP
The short-form description provided by
X.I describe_rng
can be used, at a later date, to reconstruct the state of the RNG.
This can be a slow process, since it involves seeding a virgin copy
of the random() generator, then calling random() a number of times.
The relevant call in the mrandom package is
X.I reconstruct_rng(seed, count1, count2),
where seed is the initial seed, count1 is the total number of calls to
random() since initialization, mod one billion,
and count2 is a ``div one billion'' counter.
X.LP
I hope you never have occasion to call
X.I reconstruct_rng
if count2 is greater than zero!
This is one reason to preserve occasional copies of the complete
RNG state, as recorded in the RNG state file.
If you have a state file, then
X.I restart_rng(filename)
will always give rapid results.
If the RNG state file is internally inconsistent,
restart_rng will return 0, otherwise it returns 1.
X.LP
Just before exiting a program that uses mrandom, a call to
X.I save_rng(filename)
will preserve the final state of the RNG on disk.
Note: I have found it helpful to make a
save_rng call only upon normal exit from
my experimental codes.  In case of an error exit, it is convenient
to leave the RNG state file in its initial condition.
This makes it easy to retry the experiment after fixing the bug.
X.SH AUTHOR
Clark Thomborson, [e-mail]
X.SH DIAGNOSTICS
If the error-checking code in
X.I init_rng, restart_rng,
and
X.I save_rng
discovers a problem, an error message is printed on the stderr stream.
Also, if
X.I reconstruct_rng
is called with count2 > 0, a warning message is printed on stderr.
X.LP
A demonstration routine, named
X.I mrtest.c,
is included with this software distribution to help get you started with the
X.I mrandom
package.
X.SH "SEE ALSO"
random(3), rand(3C), drand48(3)
X.SH BUGS
A little slower than
X.I random().
X.LP
As far as I know,
it is a subject of current research to decide if the pseudorandom sequence of
X.I mrandom()
is superior to that of a high-precision multiplicative generator such as
X.I drand48().
A low-precision multiplicative generator, such as
X.I rand(3C),
will cause problems in many applications, since non-random behavior
is noticeable in the least-significant twenty bits of its output.
X.LP
Currently, the best theoretical guarantees for pseudorandomness come from the
theory of universal hash functions.
The idea is to use a high-precision multiplicative generator with
random coefficients.
Thus one might consider using the mrandom package to generate random
coefficients, modulo large primes, for the drand48 pseudorandom number
generator.
X.LP
By examination of the object code for random(), I deduce that
it uses the following algorithm in its default (31-word table)
configuration.
X
X.ID
long random()
X{
X   static int j=4, k=1;
X   long r;
X   rngstate[j] += rngstate[k];
X   r = (rngstate[j]>>1) & 0x7fffffff;
X   j++;
X   if (j>31) j=1;
X   k++;
X   if (k>31) k=1;
X   return(r);
X
X}
X.DE
X
The initial values for rngstate are obtained from a linear congruential
generator.  Any defects of this scheme will be inherited by mrandom.
END_OF_FILE
if test 7553 -ne `wc -c <'mrandom.3'`; then
    echo shar: \"'mrandom.3'\" unpacked with wrong size!
fi
# end of 'mrandom.3'
fi
if test -f 'mrandom.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'mrandom.c'\"
else
echo shar: Extracting \"'mrandom.c'\" \(10569 characters\)
sed "s/^X//" >'mrandom.c' <<'END_OF_FILE'
X/*
X *        mrandom.c
X *
X *    Wrapper for random(), allowing safe restart.
X *
X *    Original Implementation:
X *        Clark Thomborson, September 1991
X *
X *    This material is based upon work supported by the National
X *    Science Foundation under grant number MIP-9023238.  The
X *    Government has certain rights in this material.
X *
X *    Any opinions, findings, and conclusions or recommendations
X *    expressed in this material are those of the author and do
X *    not necessarily reflect the view of the National Science
X *    Foundation.
X *
X *    This code is neither copyrighted nor patented.
X */
X
X#include <values.h> /* we need MAXLONG */
X#include <stdio.h>  /* we need FILE */
X#include "mrandom.h"
X
X#define RNGSTATELENGTH    32 /* legal values are 2, 8, 16, 32, 64 */
X#define RNGSTATE0    3  /* rngstate[0] for table length 32 */
X#define BILLION    1000000000 /* convenient modulus for a 32-bit counter */
X
X/* Data in our RNG statefile */
long    RNGseed;    /* the seed originally used to initialize rngstate */
long    RNGcount1;  /* mod-BILLION counter of random() calls */
long    RNGcount2;  /* div-BILLION counter */
long    rngstate[RNGSTATELENGTH]; /* an array passed to random() */
long    RNGnextval; /* random()'s next output */
X
X/* Format of our RNG statefile */
X#define RNGfileLINE1    "Initial seed = %d\n"
X#define RNGfileLINE2    \
X    "Number of mrandom() calls after seeding = %ld billion + %ld\n"
X#define RNGfileLINE3    "RNG state table =\n"
X#define RNGfileLINEn    "   %08lx %08lx %08lx %08lx\n"
X#define RNGfileLINE12    "Next value in this pseudorandom sequence = %08lx\n"
X
X/* Write a RNG state identifier into the user-supplied string rngid,
X * which must be of length at least RNGIDSTRLEN.  If the user has not
X * initialized the rng with init_rng(), restart_rng(), or reconstruct_rng(),
X * abort with an error message to stderr.  Otherwise return the value of rngid.
X */
extern char *describe_rng (rngid)
char rngid[RNGIDSTRLEN];
X{
X  if (rngstate[0] != RNGSTATE0) {
X    fprintf(stderr, "RNG has not been initialized!\n");
X    fflush(stderr);
X    exit(1);
X  }
X  sprintf(rngid, "RNG state identifier is (%ld, %ld, %ld)\n",
X         RNGseed, RNGcount1, RNGcount2);
X  return(rngid);
X}
X
X/* Create a random number statefile initialized with the given seed.
X * Return 1 if file is successfully created, 0 otherwise.
X */
extern int init_rng (seed, filename)
int seed;
char *filename;
X{
X  reconstruct_rng(seed,0,0);
X  return(save_rng(filename));
X} /* end init_rng */
X
X/* Rebuild a random() state by reseeding the generator, then making
X * (count2*1e9 + count1) calls to mrandom().  Useful for error-checking
X * and error-recovery routines, although it is very slow for large counts.
X */
extern void reconstruct_rng (seed, count1, count2)
long seed, count1, count2;
X{
X  double a;
X  initstate(seed, rngstate, RNGSTATELENGTH*4);
X  RNGseed = seed; /* keep a record of the seed */
X  RNGcount1 = 0;  /* mod-billion counter */
X  RNGcount2 = 0;  /* div-billion counter */
X  /* watch out for infinite loops */
X  if (count1 >= 0 && count1 < BILLION && count2 >= 0) {
X    if (count2 != 0) {
X      fprintf(stderr, "Warning: this reconstruction will take a LONG time!\n");
X      fflush(stderr);
X    }
X    for ( ; (RNGcount1 != count1) || (RNGcount2 != count2) ; ) {
X      a = frandom();
X    }
X  }
X} /* end reconstruct_rng */
X
X/* Return the next random() output without disturbing the RNG state.
X   This procedure will only work properly if random's state is in rngstate[].
X */
int nextval()
X{
X  long state[RNGSTATELENGTH];
X  long i, r, retval;
X  for (i=0; i<RNGSTATELENGTH; i++) {
X    state[i] = rngstate[i];
X  }
X  retval = random(); /* next value in this pseudorandom sequence */
X  for (i=1; i<RNGSTATELENGTH-1; i++) {
X    r = random(); /* these calls are necessary to restore random()'s state */
X  }
X  for (i=0; i<RNGSTATELENGTH; i++) {
X    rngstate[i] = state[i];
X  }
X  return(retval);
X}
X
X/* Restart a generator from a statefile.  Print a message on stderr
X * if the restart failed due to a garbled or non-existent statefile,
X * and return 0.  Otherwise return 1.
X */
extern int restart_rng(filename)
char *filename;
X{
X  FILE *fp;
X  long i, m, newstate[RNGSTATELENGTH], r; /* temps */
X  int errflag; /* initially 0, becomes 1 if a problem is discovered */
X
X  /* restore random()'s internal variables to their initial state */
X  m = (RNGcount1%(RNGSTATELENGTH-1)) + RNGcount2*(BILLION%(RNGSTATELENGTH-1));
X  for ( ; m%(RNGSTATELENGTH-1) != 0; m++) {
X    r = random();
X  }
X
X  /* restore counter values, retrieve original RNG seed and current state */
X  fp = fopen(filename, "r");
X  if (!fp) {
X    fprintf(stderr, "There is no RNG statefile in this directory!\n");
X    fflush(stderr);
X    exit(1);
X  }
X  fscanf(fp, RNGfileLINE1, &RNGseed);
X  fscanf(fp, RNGfileLINE2, &RNGcount2, &RNGcount1);
X  fscanf(fp, RNGfileLINE3);
X  for (i=0; i<RNGSTATELENGTH; i++) {
X    fscanf(fp, "%lx", &newstate[i]);
X  }
X  fscanf(fp, "\n");
X  fscanf(fp, RNGfileLINE12, &RNGnextval);
X  fclose(fp);
X
X  /* error check: the first word of a 32-word rngstate is always 3 */
X  if (newstate[0] != RNGSTATE0) {
X    errflag = 1;
X  } else {
X    errflag = 0;
X  }
X
X  /* tell random() we have a 32-word rngstate */
X  rngstate[0] = RNGSTATE0;
X  setstate(rngstate);
X
X  /* If reconstruction will be rapid, do it as an additional error check. */
X  if (RNGcount1 < 50 && RNGcount2 == 0) {
X    reconstruct_rng(RNGseed, RNGcount1, RNGcount2);
X    /* see if we got to the same state */
X    for (i=0; i<RNGSTATELENGTH; i++) {
X      if (newstate[i] != rngstate[i]) {
X        errflag = 1;
X      }
X    }
X  } else {
X    /* quickly modify random()'s internal state to line up with RNGcount info */
X    m = (RNGcount1%(RNGSTATELENGTH-1)) + RNGcount2*(BILLION%(RNGSTATELENGTH-1));
X    for (i=0; i < m%(RNGSTATELENGTH-1); i++) {
X      r = random();
X    }
X  }
X
X  /* copy in the new state */
X  for (i=0; i<RNGSTATELENGTH; i++) {
X    rngstate[i] = newstate[i];
X  }
X
X  /* Check nextval() operation, and verify RNGnextval */
X  if (nextval() != nextval() || RNGnextval != nextval()) {
X    errflag = 1;
X  }
X
X  if (errflag) {
X    fprintf(stderr,
X    "Warning: RNG statefile is inconsistent.  Did you edit it?\n");
X    fprintf(stderr,
X    "If not, check your program to make sure you:\n");
X    fprintf(stderr,
X    "   1. use frandom() or mrandom(m), not random();\n");
X    fprintf(stderr,
X    "   2. use init_rng(seed, filename), not srandom(seed);\n");
X    fprintf(stderr,
X    "   3. use restart_rng(filename), not setstate(state); and\n");
X    fprintf(stderr,
X    "   4. don't overwrite the private storage of mrandom.\n");
X    fflush(stderr);
X    return(0);
X  } else {
X    return(1); /* everything looks ok */
X  }
X} /* end restart_rng */
X
X/* Generate a uniformly-distributed number a, 0.0 <= a < 1.0, using
X * the 4.3bsd additive congruential generator random().
X *
X * This routine keeps track of the number of calls to random().
X *
X * From inspection of the object code for random() in a Sun-3 release,
X * I deduce the following.  A call to initstate() or setstate()
X *   1. defines the location of the state array (rngstate in the code below);
X *   2. defines the length of the state array (32 words in the code below;
X *    other possibilities are 2, 8, 16, or 64 words);
X *   3. initializes several static variables to point at the state array
X *    (the indices j and k in the code below);
X *   4. defines the randomization algorithm (a linear congruential
X *    generator is used if the state array has length 2).
X * Additionally, initstate(seed) fills the state array with the output
X * of a linear congruential generator that is seeded with the given seed.
X *
X * Here is a disassembled and decompiled listing of random(), as
X * it would operate when initialized with my init_rng routine:
X *
X * long random()
X * {
X *   static int j=4, k=1;
X *   long r;
X *   rngstate[j] += rngstate[k];
X *   r = (rngstate[j]>>1) & 0x7fffffff;
X *   j++;
X *   if (j>31) j=1;
X *   k++;
X *   if (k>31) k=1;
X *   return(r);
X * }
X *
X * Note that j and k return to their original values after every 31
X * calls to random().  This property allows me to write a restart_rng()
X * routine even though I can't read the values of j and k.
X *
X * Also note that this is a very poor RNG algorithm if it is only
X * called once between successive setstate() calls.  In this pathological
X * case, each output of random() will differ from the prior random() output
X * by rngstate[1]/2 or rngstate[1]/2+1, mod 2^31-1.  Even in non-pathological
X * calling sequences the long-period property of an additive congruential
X * generator is only guaranteed if j and k always move through the array
X * in the manner shown.  Thus the need for counting random() calls
X * (or for writing your own time-optimized additive congruential generator).
X *
X * Finally, note that integer overflow is a frequent occurrence in
X * rngstate[k] += rngstate[j].  Thus random() could conceivably
X * behave differently on different machines.
X */
extern double frandom ()
X{
X  RNGcount1++; /* mod-billion counter */
X  if (RNGcount1 == BILLION) {
X    RNGcount1 = 0;
X    RNGcount2++; /* an overflow is unlikely in my lifetime */
X  }
X  return( (double) random()/MAXLONG );
X} /* end frandom */
X
X/* Generate a random integer, uniformly distributed in the range 0..m-1.
X * We use the most-significant bits of the result of a random() call,
X * although this may only marginally improve the quality of our output.
X * You may wish to rewrite this routine if you don't have hardware
X * floating point.  If you do, be sure to include the counter-bumping
X * code from frandom().
X */
extern long mrandom (m)
long m;
X{
X  return( (int) (frandom()*m) );
X} /* end mrandom */
X
X/* Save the RNG state to a statefile, after calling random() enough
X * times to reset its internal state variables to their initial values.
X * Check to be sure the RNG can be restarted by calling restart_rng().
X * Return 0 if a problem is detected, printing an error message on stderr.
X * Otherwise return 1.
X */
extern int save_rng(filename)
char *filename;
X{
X  FILE *fp;
X  long i;
X
X  /* write the statefile */
X  fp = fopen(filename, "w");
X  if (!fp) {
X    fprintf(stderr, "Trouble opening RNG statefile %s for writing.", filename);
X    fflush(stderr);
X    return(0);
X  }
X  fprintf(fp, RNGfileLINE1, RNGseed);
X  fprintf(fp, RNGfileLINE2, RNGcount2, RNGcount1);
X  fprintf(fp, RNGfileLINE3);
X  for (i=0; i<RNGSTATELENGTH; i+=4) {
X    fprintf(fp, RNGfileLINEn,
X        rngstate[i], rngstate[i+1], rngstate[i+2], rngstate[i+3]);
X  }
X  fprintf(fp, RNGfileLINE12, nextval());
X  fclose(fp);
X
X  /* Return after checking that state was saved correctly. */
X  return( restart_rng(filename) );
X
X} /* end save_rng */
X
X/* end mrandom.c */
END_OF_FILE
if test 10569 -ne `wc -c <'mrandom.c'`; then
    echo shar: \"'mrandom.c'\" unpacked with wrong size!
fi
# end of 'mrandom.c'
fi
if test -f 'mrandom.h' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'mrandom.h'\"
else
echo shar: Extracting \"'mrandom.h'\" \(2379 characters\)
sed "s/^X//" >'mrandom.h' <<'END_OF_FILE'
X/*
X *        mrandom.h
X *
X *    A wrapper for random(), to allow safe restart.
X *
X *    Original Implementation:
X *        Clark Thomborson, September 1991
X *
X *    This material is based upon work supported by the National
X *    Science Foundation under grant number MIP-9023238.  The
X *    Government has certain rights in this material.
X *
X *    Any opinions, findings, and conclusions or recommendations
X *    expressed in this material are those of the author and do
X *    not necessarily reflect the view of the National Science
X *    Foundation.
X *
X *    This code is neither copyrighted nor patented.
X */
X
X#define RNGIDSTRLEN 80 /* max length of string written by describe_rng() */
X
char *describe_rng (/* char rngid[RNGIDSTRLEN] */);
X/* Write a RNG state identifier into the user-supplied string rngid,
X * which must be of length at least RNGIDSTRLEN.  If the user has not
X * initialized the rng with init_rng(), restart_rng(), or reconstruct_rng(),
X * abort with an error message to stderr.  Otherwise return the value of rngid.
X */
X
int init_rng (/* long seed; char *filename; */);
X/* Create a random number statefile initialized with the given seed.
X * Return 1 if file is successfully created, 0 otherwise.
X */
X
int restart_rng (/* char *filename; */);
X/* Restart a generator from a statefile.  Print a message on stderr
X * if the restart failed due to a garbled or non-existent statefile,
X * and return 1.  Otherwise return 0.
X */
X
double frandom ();
X/* Generate a uniformly-distributed number a, 0.0 <= a < 1.0, using
X * the 4.3bsd additive congruential generator random().
X */
X
long mrandom (/* long m; */);
X/* Generate a random integer, uniformly distributed in the range 0..m-1.
X * We use the most-significant bits of the result of a random() call.
X */
X
int save_rng(/* char *filename; */);
X/* Save the RNG state to a statefile, after calling random() enough
X * times to reset its internal state variables to their initial values.
X * Check to be sure the RNG can be restarted by calling restart_rng().
X * Return 0 if a problem is detected, printing an error message on stderr.
X * Otherwise return 1.
X */
X
void reconstruct_rng (/* long seed, count1, count2; */);
X/* Rebuild a random() state by reseeding the generator, then making
X * (count2*1e9 + count1) calls to random().  Useful for error-checking
X * and error-recovery routines, although it is very slow for large counts.
X */
X
X/* end mrandom.h */
END_OF_FILE
if test 2379 -ne `wc -c <'mrandom.h'`; then
    echo shar: \"'mrandom.h'\" unpacked with wrong size!
fi
# end of 'mrandom.h'
fi
if test -f 'mrtest.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'mrtest.c'\"
else
echo shar: Extracting \"'mrtest.c'\" \(2684 characters\)
sed "s/^X//" >'mrtest.c' <<'END_OF_FILE'
X/*
X *        mrtest.c
X *
X *    Test routine for mrandom.c
X *
X *    Original Implementation:
X *        Clark Thomborson, September 1991
X *
X *    This material is based upon work supported by the National
X *    Science Foundation under grant number MIP-9023238.  The
X *    Government has certain rights in this material.
X *
X *    Any opinions, findings, and conclusions or recommendations
X *    expressed in this material are those of the author and do
X *    not necessarily reflect the view of the National Science
X *    Foundation.
X *
X *    This code is neither copyrighted nor patented.
X */
X
X# include <sys/file.h> /* we use access() */
X# include "mrandom.h"
X
X# define RNGFILENAME "RNGstatefile" /* where the RNG state is stored */
X
X/* calling conventions */
void argerr(progname)
char *progname;
X{
X  printf("Usage: %s nn -snnnnnnnn -mnn R\n", progname);
X  printf("  nn sets number of random generates\n");
X  printf("  -snnnnnnnn seeds the RNG with the given value\n");
X  printf("  -mnnnn adjusts the range of the RNG\n");
X  exit(1);
X} /* end argerr */
X
X/* A simple test-driver */
int main(argc,argv)
int argc; char *argv[];
X{
X  /* command-line arguments, with default values */
X  int reseeding=0; /* nonzero if RNG state will be re-initialized */
X  int seed=1573122697; /* new seed for RNG */
X  int n=10; /* number of randoms to generate */
X  int m=100; /* desired range of random outputs: 0..m-1 */
X
X  int i; /* temp */
X  char rngdesc[RNGIDSTRLEN]; /* for a describe_rng() call */
X
X  if(argc > 1) {
X    for(i=1;i<argc;i++) {
X      if(argv[i][0] >= '0' && argv[i][0] <= '9') {
X    n = atoi(&(argv[i][0]));
X      } else if(argv[i][0] == '-') {
X    switch(argv[i][1]) {
X    case 's': /* new seed for rng */
X      seed = atoi(&(argv[i][2]));
X      reseeding = 1;
X      break;
X    case 'm': /* adjust range of rng */
X      m = atoi(&(argv[i][2]));
X      break;
X    default:
X      argerr((char*) argv[0]);
X    }
X      } else {
X    argerr((char*) argv[0]);
X      }
X    }
X  }
X
X  if (!reseeding ) {
X    if (access(RNGFILENAME, R_OK)) {
X      printf("There is no RNG statefile in this directory, so ");
X      printf("I'll make one for you.\n");
X      reseeding = 1;
X    }
X  }
X
X  if (reseeding) { /* create a new statefile from scratch */
X    if (!init_rng(seed,RNGFILENAME)) {
X      exit(1);
X    }
X    printf("Initializing RNG.   ");
X  } else { /* use an existing statefile */
X    if (!restart_rng(RNGFILENAME)) {
X      exit(1);
X    }
X    printf("Restarting RNG.  ");
X  }
X  printf(describe_rng(rngdesc));
X
X  printf("Here are %d random values in the range 0 to %d:\n", n, m-1);
X  for (i=0; i<n; i++) {
X    printf("  %d\n", mrandom(m));
X  }
X
X  /* terminate job */
X  printf("Final ");
X  printf(describe_rng(rngdesc));
X
X  return( save_rng(RNGFILENAME) );
X
X} /* end main */
END_OF_FILE
if test 2684 -ne `wc -c <'mrtest.c'`; then
    echo shar: \"'mrtest.c'\" unpacked with wrong size!
fi
# end of 'mrtest.c'
fi
echo shar: End of archive 1 \(of 1\).
cp /dev/null ark1isdone
MISSING=""
for I in 1 ; do
    if test ! -f ark${I}isdone ; then
    MISSING="${MISSING} ${I}"
    fi
done
if test "${MISSING}" = "" ; then
    echo You have the archive.
    rm -f ark[1-9]isdone
else
    echo You still need to unpack the following archives:
    echo "        " ${MISSING}
fi
##  End of shell archive.
exit 0